/*

   transmul.c -- return the multiplication of M' * M

% Version:  v0.7g
% Build:    9100818
% Date:     Oct-08 2009, 6:37 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

 */

#include "mex.h"
#include "matrix.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{

    const int *mn;
    int mnout[2], c, c1, c2, m, mc1, n, nc1;
    const double *a, *ac1, *ac1c, *ac2, *ac3;
    double *p, s, v1, v2, v3, v4, v5,
            s11, s12, s13, s14, s15,
            s22, s23, s24, s25,
            s33, s34, s35, s44, s45, s55;
    
    /* argument check */
    if (nrhs != 1 || nlhs > 1)
		mexErrMsgTxt("Bad number of input/output arguments.");
    if ((mxGetNumberOfDimensions(*prhs) != 2) ||
        (mxGetClassID(*prhs) != mxDOUBLE_CLASS))
        mexErrMsgTxt("Invalid input argument.");
    
    /* get dimensions */
    mn = mxGetDimensions(*prhs);
    m = *mn++;
    *mnout = n = *mn;
    mnout[1] = n;
    
    /* get data */
    a = (const double*) mxGetData(*prhs);
    
    /* create output matrices */
    *plhs = mxCreateNumericArray(2, mnout, mxDOUBLE_CLASS, mxREAL);
    p = (double *) mxGetData(*plhs);
    
    /* allow for special cases with small n */
    switch (n) {
        case 0: break;
        case 1:
            for (c = m, s = 0.0; c > 0; --c) {
                s += *a * *a;
                ++a;
            }
            *p = s;
            break;
        case 2:
            ac1 = &a[m];
            for (c = m, s11 = s12 = s22 = 0.0; c > 0; --c) {
                v1 = *a++;
                v2 = *ac1++;
                s11 += v1 * v1;
                s12 += v1 * v2;
                s22 += v2 * v2;
            }
            *p++ = s11;
            *p++ = s12;
            *p++ = s12;
            *p = s22;
            break;
        case 3:
            ac1 = &a[m];
            ac1c = &ac1[m];
            for (c = m, s11 = s12 = s13 = s22 = s23 = s33 = 0.0; c > 0; --c) {
                v1 = *a++;
                v2 = *ac1++;
                v3 = *ac1c++;
                s11 += v1 * v1;
                s12 += v1 * v2;
                s13 += v1 * v3;
                s22 += v2 * v2;
                s23 += v2 * v3;
                s33 += v3 * v3;
            }
            *p++ = s11;
            *p++ = s12;
            *p++ = s13;
            *p++ = s12;
            *p++ = s22;
            *p++ = s23;
            *p++ = s13;
            *p++ = s23;
            *p = s33;
            break;
        case 4:
            ac1 = &a[m];
            ac1c = &ac1[m];
            ac2 = &ac1c[m];
            for (c = m, s11 = s12 = s13 = s14 = s22 = s23 = s24 = s33 = s34 = s44 = 0.0; c > 0; --c) {
                v1 = *a++;
                v2 = *ac1++;
                v3 = *ac1c++;
                v4 = *ac2++;
                s11 += v1 * v1;
                s12 += v1 * v2;
                s13 += v1 * v3;
                s14 += v1 * v4;
                s22 += v2 * v2;
                s23 += v2 * v3;
                s24 += v2 * v4;
                s33 += v3 * v3;
                s34 += v3 * v4;
                s44 += v4 * v4;
            }
            *p++ = s11;
            *p++ = s12;
            *p++ = s13;
            *p++ = s14;
            *p++ = s12;
            *p++ = s22;
            *p++ = s23;
            *p++ = s24;
            *p++ = s13;
            *p++ = s23;
            *p++ = s33;
            *p++ = s34;
            *p++ = s14;
            *p++ = s24;
            *p++ = s34;
            *p = s44;
            break;
        case 5:
            ac1 = &a[m];
            ac1c = &ac1[m];
            ac2 = &ac1c[m];
            ac3 = &ac2[m];
            for (c = m, s11 = s12 = s13 = s14 = s15 = s22 = s23 = s24 = s25 = s33 = s34 = s35 = s44 = s45 = s55 = 0.0; c > 0; --c) {
                v1 = *a++;
                v2 = *ac1++;
                v3 = *ac1c++;
                v4 = *ac2++;
                v5 = *ac3++;
                s11 += v1 * v1;
                s12 += v1 * v2;
                s13 += v1 * v3;
                s14 += v1 * v4;
                s15 += v1 * v5;
                s22 += v2 * v2;
                s23 += v2 * v3;
                s24 += v2 * v4;
                s25 += v2 * v5;
                s33 += v3 * v3;
                s34 += v3 * v4;
                s35 += v3 * v5;
                s44 += v4 * v4;
                s45 += v4 * v5;
                s55 += v5 * v5;
            }
            *p++ = s11;
            *p++ = s12;
            *p++ = s13;
            *p++ = s14;
            *p++ = s15;
            *p++ = s12;
            *p++ = s22;
            *p++ = s23;
            *p++ = s24;
            *p++ = s25;
            *p++ = s13;
            *p++ = s23;
            *p++ = s33;
            *p++ = s34;
            *p++ = s35;
            *p++ = s14;
            *p++ = s24;
            *p++ = s34;
            *p++ = s44;
            *p++ = s45;
            *p++ = s15;
            *p++ = s25;
            *p++ = s35;
            *p++ = s45;
            *p = s55;
            break;
        default:
            /* loop over columns (twice) */
            for (c1 = 0; c1 < n; ++c1) {
                c2 = c1;
                mc1 = c1 * m;
                nc1 = c1 * n;
                ac1 = &a[mc1];
                for (c = m, s = 0.0, ac2 = ac1; c > 0; --c) {
                    s += *ac2 * *ac2;
                    ++ac2;
                }
                p[nc1 + c1] = s;
                for (++c2; c2 < n; ++c2) {
                    for (c = m, s = 0.0, ac1c = ac1; c > 0; --c)
                        s += *ac1c++ * *ac2++;
                    p[nc1 + c2] = p[c2 * n + c1] = s;
                }
            }
            break;
    }
}
